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SMALL  UNSTEADY  PERTURBATIONS  IN  TRANSONIC  FLOWS+ 

* **  *** 

K-Y.  Fung  , N.  J.  Yu  and  R.  Seebass 

University  of  Arizona 
Tucson,  Arizona  85721 


Abstract 

We  investigate  the  effects  of  very  small,  low  frequency,  perturbations 
to  steady  transonic  flows.  We  do  so  in  the  context  of  two-dimensional  flows 
described  by  the  small  perturbation  equation.  We  draw  inferences  from  an 
even  simpler  model  equation.  Our  primary  concern  is  with  the  validity  of 
linearizing  the  unsteady  perturbations  to  such  flows  and,  in  particular, 
with  the  failure  of  earlier  studies  to  account  for  the  shock  wave  motions 
that  we  know  occur.  We  provide  a method  that  allows  one  to  account  for 
shock  wave  motions  due  to  arbitrary,  but  small,  unsteady  changes  in  the 
boundary  conditions.  Consequently,  both  harmonic  and  indicial  responses 
may  be  determined.  Time-linearized  results  for  the  transonic  flow  past  an 
NACA  64A006  airfoil  experiencing  harmonic  motions  in  one  of  several  modes 
are  presented.  Selected  results  are  compared  with  those  obtained  from 
nonlinear  calculations  using  a shock-fitting  algorithm. 
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Introduction 


In  unsteady  transonic  flows,  relatively  small  periodic  changes  in  the  1 

boundary  conditions  can  lead  to  substantial  changes  in  the  magnitude  and  i 

phase  lag  of  loads  and  moments.  These  are  of  major  concern  in  the  aerody-  I 

namic  design  of  aircraft  that  operate  in  the  transonic  regime.  Reference  1 
contains  a short,  but  timely,  review  of  various  aspects  of  unsteady  tran- 
sonic flow.  Of  particular  concern  are  aeroelastic  behavior  and  flutter 

boundaries.  Here  the  unsteady  perturbations  may  be  considered  small,  and 

2 

linearization  about  a nonlinear  steady  flow,  as  suggested  by  Landahl  long 
ago,  would  seem  to  be  appropriate.  Difficulties  arise,  however,  that  detract 
from  this  procedure.  While  the  equation  is  linear,  its  coefficients  are 
variable  and  must  be  determined  by  numerical  solution  of  a nonlinear  problem 
that,  in  the  cases  of  prime  interest,  has  a discontinuous  solution;  that  is, 
there  are  embedded  shock  waves.  Also,  while  a change  of  variables  in  the 
linear  equation  provides  a scaling  of  parameters  that  is  indicative  of  the 
trade-offs  between,  e.g.,  Mach  number  and  reduced  frequency,  the  only  simili- 
tude is  the  one  basic  to  the  nonlinear  formulation. 

3 4 

Traci  et  al.  ’ have  developed  relaxation  methods  for  solving  the  result- 
ing time-linearized  equations  of  motion.  Less  complete,  but  comparable, 
studies  have  been  made  by  Weatherill  et  al.^;  these  derive  from  an  earlier 

study  by  Ehlers^.  In  both  of  these  studies  shock  motions,  which  contribute 

7 8 

substantially  to  the  time  varying  loads  and  moments,  ’ are  neglected.  Also, 

difficulties  arise  in  the  convergence  of  the  iterative  numerical  scheme. 

9 

Here  we  pursue  a different  numerical  course.  Yu  et  al.  have  developed 
a numerical  procedure  for  computing  solutions  to  the  unsteady  small  perturba- 
tion equation  for  transonic  flows,  which  treats  embedded  shock  waves  as  dis- 
continuities. This  procedure  can  be  used  to  calculate  the  basic  steady  flow 
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that  we  wish  to  subject  to  small  unsteady  perturbations.  A simplified  version 
of  this  algorithm  can  then  be  used  to  calculate  the  linearized  unsteady  per- 
turbations to  the  flow.  These  calculations  can  be  carried  out  in  conjunction 
with  an  algorithm  that  determines  the  shock  wave  motion.  The  procedure  we 

have  used  to  calculate  the  shock  wave  motion  is  a rather  obvious  one;  it  is 

2 

not  surprising  then,  that  it,  too,  was  given  in  the  monograph  by  Landahl 
(Section  10.2).  An  alternative  procedure,  related  in  some  ways  to  that  used 
here,  is  implied  by  Nixon's^  study  of  perturbations  to  steady  discontinuous 
transonic  flows. 

After  a brief  review  of  the  formulation  of  time-linearized  methods,  we 
discuss  a one-dimensional  model  equation.  This  provides  a suitable  testing 
ground  for  our  ideas  and  serves  to  illustrate  several  basic  points  we  wish  to 
stress.  Next,  we  discuss  the  two-dimensional  formulation  and  compare  results 
of  time-linearized  calculations  with  those  obtained  without  the  linearization, 
for  an  NACA  64A006  airfoil  oscillating  in  pitch. 

Formulation 

We  write  the  unsteady  small  disturbance  equation  for  low  frequency  tran- 
sonic flows  in  the  commonly  used  form 

-2KM2<J>  _ + {1  - M2  - (y  + 1)M24>  }<fr  + <p  - 0.  (1) 

» xt  ® ® x XX  yy 

The  spatial  coordinates,  the  time,  and  the  velocity  potential  in  (1)  have  been 
non-dimensionalized  by  the  chord,  the  reciprocal  of  the  angular  frequency,  and 
the  free  stream  velocity  times  the  chord,  respectively.  Other,  perhaps  more 
suitable,  e ^s  are  given  in  Reference  9.  This  equation  results  from  a 
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r 


systematic  expansion  of  the  velocity  potential  in  the  thickness  ratio  r 


2/3, 


and  applies  for  reduced  frequencies  K = 0(t  ),  where  K = ojc/U,  i.e., 


the  angular  frequency  multiplied  by  the  time  it  takes  the  flow  to  traverse 

11 


the  airfoil  chord.  Lin,  Reisner  and  Tsien  showed  that,  with  restriction 
to  small  perturbations  throughout  the  flow,  (1)  is  the  only  nonlinear  equa- 
tion that  arises.  For  moderate  frequencies  the  equation 


-K2*tt  - 2K*xt  + {1  - M2  - (y  + 1)M2[*x  + K^]}^  + « 


yy 


is  frequently  used,  with  or  without  the  $ term,  and  may  provide  results 


that  apply  at  higher  frequencies  than  those  obtained  from  (1)  or  the  linear 

3 


form  of  the  above  equation. 

2/3 


Because  K * 0(t  ),  the  boundary  condition  on  the  body  takes  the  sim- 

ple form 


<t  (x,0,t) 


t3Y(x, t)/3x  - r[Y°  + ~(YU  + Ky“)], 

X T X C 


(2) 


where  Y(x,t),  the  instantaneous  body  shape,  has  been  decomposed  into  a 


u u 

steady  part,  Y , and  an  unsteady  part,  Y . The  last  term,  KYc>  is  dropped 


except  when  Y^  is  small  or  zero  because  K = 0(t2^).  Here  6 is  the  amp- 


litude of  the  unsteady  oscillation.  Far  from  the  body  we  require  that  the 
derivatives  of  $ vanish.  In  this  approximation  the  pressure  coefficient, 


defined  so  that  it  vanishes  at  sonic  conditions,  takes  the  form 


/ 


C - -2  < 
P ' 


M2  - 1 


(Y  + DM 


2 + *x 


(3) 


In  the  small  disturbance  approximation,  the  Kutta  condition  is  imposed  by  re- 


quiring that  Cp  be  continuous  at  y » 0 for  x > 1/2. 
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Any  shock  wave  that  exists  in  the  flow  field  must  satisfy  the  jump  rela- 
tion derived  from  the  conservative  form  of  the  governing  equation  (1),  namely, 

-2KMf|[*xl|2(^)  - U - M2  - (y  + DM2; JU J2  + l4»  I2  « 0 (4) 

together  with  the  condition  derived  from  the  assumption  of  irrotationality , 

““1+xl/lO-  (5) 

s y 

Here  refers  to  the  mean  value  of  <J>  evaluated  on  each  side  of  the  dis- 

^ X 

continuity,  and  indicates  the  jump  in  <j>x  across  the  discontinuity; 

the  subscript  "s"  denotes  the  quantity  evaluated  at  the  shock  surface. 

Time-Lineari2ed  Equations 

We  now  assume  that  the  unsteady  disturbances,  characterized  by  6,  are 
small  enough  that  we  may  write 


4>(x,y,t)  “ <p°  (x,y)  + 5i(;(x,y,t)  + o(<5) 


(6) 


and  neglect  higher-order  terms  in  5.  The  restriction  imposed  on  6 for  this 

to  be  true  will  depend  on  the  other  parameters  of  the  problem,  viz., 

2 2 2/1 

k = (1  - + DM^t]  and  K.  This  gives 


U - 


00 


1 

- 2±x 


(7) 


and 


♦J(x,0)  « tY° ' (x) , 


The  solution  to  (7)  must  satisfy  the  steady  version  of  the  shock  relations 
(4)  and  (5).  The  shock  relations  for  (8)  are  discussed  later. 


As  mentioned  above,  a shock-fitting  scheme  that  approximates  the  shock 

9 

waves  as  discontinuities  normal  to  the  free  stream  has  been  developed  with 

an  alternating-direction  implicit  scheme  (i.e.,  ADI)  to  compute  the  solution 

12 

to  (7).  Comparison  of  these  results  with  the  results  obtained  using  an 
exact  shock-fitting  algorithm  and  line  relaxation  indicate  that  they  should 
suffice  for  most  studies.  At  the  very  least  they  should  prove  adequate  for 
the  time-linearized  studies  of  interest  here,  as  only  small  shock  excursions 
can  be  allowed. 

We  assume,  then,  that  we  have  the  numerical  values  for  <j>^  required  in 
(8).  These  data  will  be  discontinuous  across  some  vertical  line,  the  shock 
wave,  x = x*,  0 <_  |y|  <_  y*.  We  then  ask,  under  what  conditions  is  (8) 

valid?  And  how  do  we  account  for  shock  wave  motions  in  the  linearized  analy- 
sis? The  answers  to  these  questions  are  inferred  from  a simple  one-dimensional 
model  discussed  in  the  next  section. 

Anticipating  that  we  will  wish  to  solve  (8)  with  the  same  technique  that 
proved  successful  for  (1)  we  avoid  writing 

tKx,y,t)  =“  Re{iKx,y)eia)t}.  (9) 

The  assumption  (9)  restricts  the  study  to  harmonic  linear  motions.  Because 
indicial  motions  can  be  superimposed  to  obtain  the  results  for  any  frequency  , 


■■■mb mmm  i 
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they  too  seem  important.  The  assumption  (9)  suppresses  the  time  dimension  of 

the  calculation  but  it  results  in  two  coupled  equations,  or  one  equation  for 

a complex-valued  ip,  which  may  be  solved  by  line  relaxation.  Our  experience 

with  unsteady  ADI  techniques  has  been  that  they  are  at  least  as  effective  as 

line  relaxation  for  problems  of  this  type,  and  hence  there  is  no  numerical 

advantage  to  the  decomposition  (9).  This  conclusion  was  also  arrived  at  by 
13 

Ballhaus  et  al.  in  a related  study. 

An  appropriate  scaling  of  the  dependent  and  independent  variables  in  (7) 
and  (8)  allows  either  the  thickness  or  the  frequency  to  be  normalized  to  the 
value  1,  as  expected.  This  scaling,  in  terms  of  the  transonic  similarity 
parameter  k,  the  amplitude  of  the  unsteady  motion  6,  its  frequency  gj, 
and  the  body's  basic  thickness  t,  leads  to 


iKx,y;<;6;u;l) 


uip 


x 


y k ..  1 

’ 1/2’  l/2;i;  1/2 

U)  0)  OJ 


/ 


where  x and  y are  suitably  scaled  replacements  for  the  x and  y coordi- 
nates. This  result  can  be  used  to  check  trends  noted  in  the  numerical  results. 


One-Dimensional  Model 

To  answer  the  questions  raised  above,  we  study  a simple  unsteady  one- 
dimensional analog  of  (2).  We  consider  a one-dimensional  unsteady  equation 
that  models  the  important  features  of  (2),  and  ascertain  how  a simple  steady 
solution  is  modified  by  small  unsteady  perturbations.  We  consider,  then, 


2*xt  + (1  - 


> H 

X XX 


2<f> 


xt 


> )2i 

X X 


(10) 


subject  to  <J>  (0 , t ) = f^Ct),  d>x  (0,  t ) - f 2 ( t)  and  either  $(l,t)  or  4>x  (1 , t ) * 
f(t).  There  are  restrictions  on  f which,  for  brevity,  we  do  not  list.  Our 
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study  could  be  generalized  by  replacing  1 — d>  in  (10)  by  (A(x)  - $ ) 

x x 

where  A(x)  is  a continuous  function  of  x,  but  little  added  insight  is 
gained. 

To  simplify  matters  further,  consider  the  especially  simple  subcase 
^2  * 0»  ^2  = f = 3 + 5p(t).  When  6=0,  we  have  the  steady  solution 

-x,  0 <_  x <■  3/4 

-3(1  - x) , 3/4  £ x <_  1. 

This  satisfies  (10)  and  the  jump  condition  that  one  derives  from  it,  viz., 

dx 

I<M  = 0 on  2 = 1 - 4»x>  (11) 

Now  a general  solution  of  (10),  in  terms  of  <f>x>  is 
<f>x  = arbitrary  function  of 

This  can  be  verified  by  substitution.  With,  say  <j>x(l,t)  = 3 + 6p(t)  we 

have,  for  x > x , 
s 

<{>(x,t)  ■ 3(x  - x ) + 6 f p t + ~2(1  ~ ? — dx  - h(t)  (12) 

Xs  l 1“*x^t>/ 

where  the  choice  h(t)  = x (t)  insures  that  I4>1  - 0 at  x * x because 

s s 

$ * for  x < x . 

s 

The  shock  motion  must  be  determined  by  the  direct  integration  of  (11): 
dx  . 

2 df  = 1 " *x  = ' 2 


(13) 
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Now,  for  comparison,  we  determine  the  results  that  are  obtained  by  time 
linearization;  i.e.,  we  write 


4>(x,t;<5)  " 4>°(x)  + Sip  (x,  t)  + o(6) 


and  solve  the  linear  equation  for  ip  that  results  by  dropping  higher-order 
terms  in  5.  That  is,  we  solve 


2ip  + (1  - <p°)ip  - <p°  ip  =0 
Xt  X XX  XX  X 


subject  to  i/>x(l,t)  = p(t). 

We  now  linearize  the  first  of  equations  (11)  as  follows: 


I<t>(x,t)J  * 0 * iU(x°,t)  + 4)  (x°,t)(x  (t)  - x°)  + ...  1 

s X s s s 


>°(x°)  + SiKx°,t)  + <c°(x  (t)  - x°: 

S S X s s 


Thus  we  conclude  that 


<5  (*(x",t)]  = -IU°(x0)J  (x  - x°) 

S X s s s 


From  the  second  equation  of  (11),  with  x (t)  * x°  + <$x(t)  we  find 

s s 


dt  2 ^x  4 x^ 


where  ( )^  refers  to  the  value  behind  the  shock.  Thus  we  may  replace  (11) 


■I 


-10- 


!♦(*•,  t)l  = -I<(x»)]x 

S X s 

or 


and 


i|/  (x®,t)  = -4x(t), 

D S 


(17a) 

(17b) 


Example 

Consider  now,  for  example, 

p(t)  = sin  iot; 

it  is  easy  to  show  that  a general  solution  to  (15)  is 


(18) 

(19) 

Had  we  solved  (15)  with  ip  = 0 for  x < x (t)  and  determined  the 

s 

exact  shock  motion  from  (11)  with  1 — «(>  — —2  +0(6)  used  in  (12),  we 

would  find  that  behind  the  shock 


<Kx,t)  = [cos  “(1  - x - t)  - h(t)  ] . 

The  function  h(t)  follows  from  (17a, b)  and  assuming,  e.g.,  that 
i|>(3/4, 0)  =*  0.  Thus 

i/>(x,t)  = - ^ [cos  io(l  - x - t)  - cos  (<o/4)  ] 

and 

x(t)  - ^ [cos  <o(-£  - t)  - cos  (to/ 4)  ] . 


$ (x , t ) = <(>0  + 3-4x  - - [cos  U)(l  - x - t)  - cos  oo(l  - x - t ) ] , 

3 (jJ  s 


(20) 
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where  the  shock  motion  is  given  explicitly  by 


2 _x  I tan  | (t  + c)  - [tan  j (t  - 1)  - -|] 

x » — tan  ( — 

s w , , w 0)  , x <5, 

1 + tan  — (t  - 1 ) [ tan  — (c  + t)  + — ] 
2 2 4 


with  c=-2(tan  [tan  (u)/8)  + 6/4]  )/o).  For  o>  >>  6,  which  is  required  for 


small  shock  motions,  (21)  simplifies  to 


3 (S  1 

x = — + - — [cos  u>(y-  - t)  - cos  (w/4)]. 
s 4 4o)  4 


Thus  these  results  are  in  agreement  with  (18)  and  (19)  to  lowest  order  in  <5. 

The  time-linearized  results  (18)  and  (19)  are  now  compared  with  the  exact  results. 

The  nonlinear  result  for  x > x , given  by  (12)  and  (13),  is 

s 


+ 3 - 4x  (t)  + 5 sin  u t + 


2(1  - x) 

1 - 4>x(x,t) 


where 


dx  , / 1 - x 

s _ • 5 ^ . s 

-j—  = x » — — sin  <d  t + 

dt  s 4 . -• 

1 - 2x 


The  results  (22)  and  (23)  are  consistent  with  the  time- linearized  re- 
sults (18)  and  (19)  to  0(6),  except  for  a slow  secular  drift  in  the  shock 
2 

position  of  0(5  t)  that  occurs  in  (23)  but  not  in  (19).  This  is  an  arti- 
fact of  our  one-dimensional  model;  even  if  it  were  not,  it  would  not  invali- 
date the  use  of  the  linear  results  for  flutter  studies  where  5 is  small  and 
structural  damping  determines  the  time  scale  of  interest. 

The  time-linearized  result  for  <J,  given  by  (14)  and  (18),  is  shown  in  Figure  1 
for  5 = 0.1  with  w = 0.5.  The  shock  motion  (19)  is  shown  in  Figure  2, 
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again  for  6 * 0.1  and  with  u K 0.5,  1.0  and  2.0.  As  is  clear  from  (19), 

the  shock  excursions  are  0(6/w).  Consequently,  they  cannot  be  neglected  in 

the  calculation  unless  6/w  is  much  smaller  than  one.  For  the  times  shown, 

there  is  no  significant  difference  between  the  results  (19)  and  (23)  due  to 
2 

the  0(6  t)  term. 

The  main  conclusions  we  derive  from  this  study  are  that  it  is  essential 
to  consider  shock  motion  in  computing  time-linearized  solutions  if  we  are 
to  determine  the  effects  of  small  unsteady  perturbations  correctly  to  lowest 
order,  and  that  shock  excursions  increase  as  the  frequency  is  decreased. 
Additionally,  this  motion  can  be  computed  in  a straightforward  manner. 

Two-Dimensional  Time-Linearized  Analysis 

The  results  from  our  simple  model  show  that  the  time-linearized  results 
must  be  corrected  for  shock  motions  if  they  are  to  be  consistently  correct  to 
lowest  order.  This  can  be  accomplished  by  calculating  the  shock  motion  in 
conjunction  with  the  time-linearized  solution.  Here  we  follow  an  analogous 
procedure  and  calculate  the  change  with  time  of  the  values  of  the  perturbed 
potential  behind  the  shock  required  by  the  linearized  shock  jump  relations. 
Thus,  we  wish  to  solve  (8) 

-2KMJJ*  + {[1  - M*  - (Y  + 1)mV]*}  + * = 0 

oo  xt  °°  00  x x x yy 

with  (8) 

I{iy(x,0,t)  = Y^(x,t)  * “ \ 1 x £ \ 

subject  to  the  far-field  boundary  and  Kutta  conditions.  As  we  noted,  the 
steady  result  for  4>°  can  be  calculated  adequately  for  most  small  disturbance 


flows  using  normal  shock  fitting  as  described  in  Reference  9.  Under  the 
assumption  that  the  shock  wave  is  normal,  the  shock  jump  conditions  (4)  and 
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(5)  can  be  replaced  by  requiring 


dx 


M2-  1 


(y  + 1)M, 


2 Tx 


(24) 


For  steady  flow  xg  =•  0 and  in  (24)  {•••}  = 0.  We  express  the  shock  posi- 

tion as 


x * x°  + 5x(t) 
s s 


and  conclude  that  the  shock  motion  is  governed  by 


& = x.+ 1 ; 

dt  2K  V 


As  discussed  in  Reference  9,  <p ^ is  evaluated  at  y = 0.  On  the  shock 


I d>I  - + 614-1; 


(25) 


linearizing  the  expression  (25)  for  the  velocity  potential  about  the  steady 
shock  position  we  find 


$(x  ,y,t)  - <P(x°, y,t)  + 4>  (x°,y,t)6dx 

s S X s 

- $°(x8  y)  + r<x°,y)6dx  + «*(x°,y,t)  + 0(52). 

s X s s 

Because  we  have  treated  the  shock  as  a normal  one,  y appears  here  simply  as 

a paramet^i. . Now  1>(*  , t,y)3  and  I$8(x°,y)J  are  both  zero;  consequently 

s s 


we  have 
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I ip  (x°  ,y , t)]| 
s 


(y  + 1) 
2K 


^(x;,y)l 


,t_  - . 

'If  (x°,0,t)dt 
0 x s 


(26) 


which  must  be  integrated  in  time  in  conjunction  with  the  solution  to  (8). 

Equation  (8)  is  now  solved  numerically  in  time  and  space  in  conjunction 
with  (26),  which  is  used  to  update  the  values  of  ip  behind  the  shock.  We 
start  with  a steady  solution  and  initiate  a body  motion,  such  as  the  harmonic 
oscillation  of  a flap.  The  calculations  proceed  in  time  until  they  are 
judged  to  be  periodic.  Note  that  indicial  as  well  as  harmonic  motions  may 
be  considered  because  we  have  not  utilized  the  usual  harmonic  decomposition 
(9). 


1 

Numerical  Procedure 

The  numerical  procedure  used  here  derives  from  that  developed  for  the 
nonlinear  equation  (1).  The  main  simplification  occurs  in  the  shock  jump  con- 
ditions. In  order  to  minimize  the  far-field  boundary  effects  on  the  results 
14 

which,  as  Magnus  has  noted,  can  be  significant,  we  again  use  coordinate 
9 

stretching  in  the  form: 


C = ±[1  - exp  (— a^ | x | ) ] for  x < 0, 
n = ±[1  - exp  (-a  |y|)]  for  y < 0, 


where  a^  and  a^  are  constants  that  determine  the  mesh  distribution.  This 
stretching  transforms  the  infinite  physical  domain  into  a computational  domain 
bounded  by  |c|  <_  1 and  | n | <_  1.  The  mesh  distribution  is  concentrated  on 
the  airfoil.  In  these  coordinates  (8)  becomes 


Wt  + A2{f(5’n)(1  " lCl>Vc  + A3{(1  “ 


ultk 


(27) 
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where 


-2K&i 

a2(l  - I n [ ) ’ 


A. 


a,(l  - | n| ) ’ 


1 - 5 


and 


f (5»n)  - 1 - M*  - ax(Y  + 1)M*(1  - |C|H° 


is  known  in  discrete  form  from  the  steady  numerical  solution.  This  function 
is  discontinuous  at  C = for  0 < n < n*.  First-order  backward  time  and 
spatial  differences  are  used  for  the  first  term.  Centered  or  first-order 
backward  differences  are  used  for  the  second  term  if  f is  less  than  or 
greater  than  zero,  respectively;  f(£,n)  is  known  in  advance  with  the  deri- 
vative automatically  evaluated  correctly.  Centered  differences  are 

used  for  the  third  term  and  denoted  by  5 . 

n 

The  solution  is  computed  using  an  alternating-direction  implicit  pro- 
cedure first  applied  to  transonic  flow  problems  by  Ballhaus  and  Steger^  and 
by  Beam  and  Warming^,  and  subsequently  further  refined  by  Ballhaus  and 
Goorjian^.  The  solution  is  advanced  in  time  from  its  initial  steady  state 
to  subsequent  time  levels  with  the  following  two-step  procedure. 

New  values  of  ij>,  denoted  by  ij>+,  are  calculated  along  n “ constant 
lines  using 


i>+  - 

A1  * + A2{f(5,n)(l  - + A36n{(1  " M >♦”*  * °* 


This  is  coupled  with  the  computation  of  new  values  of  behind  the  shock 

obtained  by  using  (26).  With  the  shock  located  at  ?°  such  that 

s 


« < c!  < 5 


g+^»  we  can  express  the  values  of  i>  ahead  of  and  behind  the 
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shock  in  a Taylor  series,  finally  arriving  at  the  result 


I^l  + = -c(n)Att^  + IUln 


(28) 


where 


c(n) 


and  <l>  is  evaluated,  following  (26),  at  n = 0.  One-half  the  change  in 
ip  across  the  shock  is  accounted  for  in  this  step,  effectively  using  the 
trapezoidal  rule  in  the  time  integration  (26);  hence  C is  half  the  value 
implied  by  (26). 

With  the  values  of  i|j+  determined,  the  new  values  of  ip  at  the  sub- 
Q+l 

sequent  time  level,  ip  , are  calculated  using 
^n+l  ^ 

Ai  £ + r V(1 ' N>‘C  - ■ 0 


in  conjunction  with  the  completion  of  the  time  integration  (26), 


I*Jn+1  - -c(n)At^+1  + I*J+. 


(29) 


~n+l  * 

Again,  ip  is  evaluated  at  n ■ 0. 

The  full  procedure  is,  effectively, 


i|»n+1  - i|>n 

Ax  ~L-  L + A2{f(C,n)(l  - |5|)^>5  + \ A36ri((l  - |ri|)0l^+1  + ♦")} 


The  results  given  by  the  authors  in  AIAA  Paper  No.  77-675  were  in  error 
because  <p.  was  allowed  to  vary  with  n;  this  is  not  consistent  with  the 
normal  shodk  approximation  that  gives  (24). 
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with  (26)  implemented  in  the  form  (28),  (29).  The  procedure  outlined  here 
effectively  corrects  the  i Ji  values  for  shock  motions  as  the  solution  pro- 
gresses. The  shock  motion  is  easily  determined  simultaneously  by  using  (26) 
and  the  expression  for  dx/dt  to  find 

Xn+1(0,t)  - -I<Kx°,0,t)]n+1/IUx(x°,0)]|.  (30) 

The  computations  then  provide  results  for  <f>x  like  those  sketched  in 
Figure  3.  This  figure  depicts  the  steady  state  result  and  the  unsteady 
changes  as  well  as  the  shock  positions  at  two  different  time  levels  where  the 
shock  is  behind  the  steady  state  position.  When  the  shocks  have  been  inserted 
in  their  known  positions  we  see  that  we  need  to  analytically  continue  data 
ahead  of  and  behind  the  shock  in  order  to  complete  the  solution.  For  shock 
excursions  that  are  o(l)  we  can  simply  extrapolate  the  steady  state  data, 
both  ahead  of  and  behind  the  shock,  to  determine  the  pressure  distribution 
on  the  body  correctly  to  lowest  order.  Larger  shock  motions  are,  of  course, 
not  admissable  in  this  theory. 


Results  and  Discussion 

Time-linearized  results  have  been  computed  for  an  NACA  64A006  airfoil 
experiencing  harmonic  pitching  and  flap  motions.  As  noted  earlier,  in  the  low 
frequency  approximation  made  here, pitching  and  plunging  motions  lead  to  the 
same  result  except  that  the  time-linearized  perturbations  are  proportional  to 
the  maximum  pitch  angle  for  the  former,  and  K times  the  maximum  amplitude 
for  the  latter.  Harmonic  motions  initiated  from  a steady  state  become  nearly 
periodic  in  three  to  ten  cycles,  with  the  changes  induced  by  flap  oscillations 
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becoming  periodic  more  rapidly  than  those  resulting  from  pitching  oscilla- 
tions. More  cycles  were  required  for  larger  reduced  frequencies  and,  to  a 
lesser  degree,  higher  Mach  numbers. 

In  order  to  confirm  the  validity  of  the  time-linearized  calculations, 
both  the  time-linearized  and  nonlinear  algorithms  were  used  to  compute  the 
response  to  a step  change  in  angle  of  attack  and  the  harmonic  response  to 
pitching  motions.  Figure  4 compares  the  nonlinear  and  time-linearized  re- 
sults for  the  normalized  circulation  and  shock  position  for  harmonic  pitching 
motions  at  * 0.88  and  K =*  0.48.  For  these  conditions  very  small  un- 
steady changes  lead  to  very  small  shock  motions  and  in  both  calcu- 
lations the  shock  wave  remains  between  grid  points.  Because  of  the 
extrapolation  procedure  used  in  the  nonlinear  shock-fitting,  the  finite  mesh 
size  used  can  introduce  errors,  albeit  small  ones,  in  the  shock's  position 
when  a grid  line  is  crossed.  We  wished  to  eliminate  these  errors  in  order  to 
use  the  nonlinear  calculations  to  judge  the  accuracy  of  time-linearized  cal- 
culations. These  results  indicate  that  for  pitching  about  mid-chord,  non- 
linear, amplitude  dependent,  behavior  occurs  for  6/t  0.1  for  K = 0.48. 

Because  the  amplitude  of  the  shock  motion  increases  with  decreasing  K,  non- 
linear effects  occur  at  smaller  values  of  5/t  at  lower  reduced  frequencies. 
Results  are  given  for  the  fifth  cycle;  note  that  the  nonlinear  results  are 
not  yet  periodic.  Figure  5 compares  the  nonlinear  and  time-linearized  pres- 
sure deviation  from  steady  state  at  six  angular  times  for  the  same  conditions. 
Good  agreement  between  the  results  is  obtained  for  <5/x  less  than  0.1. 

Time-linearized  pressure  distributions  at  six  angular  positions  for  an 
oscillating  quarter-chord  flap  with  K ■ 0.06  and  M^  * 0.875  are  shown  in 
Figure  6.  The  flap  deflection  is  downward  during  the  first  half  of  the  cycle. 
The  results  for  the  second  half  of  the  period,  for  the  symmetrical  problem 


i 
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shown  here,  are  just  the  results  shown  with  the  lower  and  upper  surface  pres- 
sures interchanged.  Thus  the  results  for  0°  are  not  given  as  they  are  just 
those  for  180°  with  the  lower  and  upper  surface  pressures  reversed.  Be- 
cause the  flap  hinge  occurs  very  close  to  the  steady  state  shock  location, 
the  pressure  singularity  due  to  the  change  in  flow  direction  at  the  hinge  is 
missed.  The  circulation  and  shock  excursion  obey  the  following  relations: 

r(t)/<5  - 9.26  sin  (t  - 59°), 

X (t)  - 12  sin  (t  - 51°). 

Note  the  substantial  phase  lag  in  the  circulation  and  the  shock's  position. 

Time- linearized  pressure  distributions  at  six  angular  positions  for  an 
oscillating  airfoil  with  K * 0.12  and  ■ 0.875  are  depicted  in  Figure 
7.  If  these  results  are  multiplied  by  K,  then  they  represent  the  pressure 
perturbations  for  a plunging  airfoil.  As  in  the  previous  case  of  an  oscil- 
lating flap,  changes  in  forces  and  moments  of  0(6/K)  occur  due  to  shock 
wave  motion.  In  this  case 

T(t)/<S  - 5.48  sin  (t  - 70°), 

X(t)  - 5.62  sin  (t  - 87°). 

Analogous  computations  have  been  carried  out  for  K ■ 0.06,  0.12,  0.24 
and  0.48.  Figure  8 depicts  the  shock  wave's  excursion  and  maximum  circulation 
as  a function  of  K The  nearly  linear  variation  of  the  shock  excursion 
substantiates  an  observation  made  in  a one-dimensional  model  where  the  shock 
wave  excursion  is  directly  proportional  to  1/K. 
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In  these  calculations  the  circulation  gives  an  immediate  evaluation  of 

the  lift  coefficient  as  a function  of  time;  the  moment  coefficient  must  be 

evaluated  by  integrating  the  moment  of  the  pressure  coefficient.  This  is 

done  by  integrating  the  moment  of  pressure  perturbations  with  the  shock  wave 

in  its  steady-state  position  and  then  correcting  these  results  for  the  moment 

due  to  the  shock  wave  motion,  assuming  that  the  shock's  strength  is  defined 

by  the  steady-state  pressure  field.  This  makes  an  error  in  the  shock  strength 

2 

of  0(6),  but  the  effect  on  the  moment  is  0(6  /K) ; because  we  have  neglected 
other  higher-order  terms  it  is  consistent  to  neglect  this  change  in  the 
strength  of  the  shock  wave. 

For  the  time-! * near! zed  results  to  be  valid  we  must  really  require 
6/tK  <<  1.  Our  numerical  results  indicate  that  for  6/xK  <_  0.2  the  unsteady 
perturbations  are  essentially  linear. 

The  time- linearized  algorithm  used  here  is  a derivative  of  that  used  for 
the  nonlinear  calculations.  Consequently,  computational  times  are  not  greatly 
reduced  from  those  required  for  the  nonlinear  calculations.  The  linearity  of 
these  computations  may  make  it  possible  to  greatly  reduce  the  computational 
effort  required.  Numerical  experience  has  shown  some  difficulties  for 
At (in  degrees) /K  _>  50.  This  is  in  agreement  with  the  consistency  requirement 
for  the  ADI  algorithm  used  here.  Both  the  domain  of  dependence  condition 
and  a local  linearized  stability  analysis  shows  the  procedure  to  be  uncondition- 
ally stable.  Each  time  step  requires  about  two  seconds  of  CPU  time  on  a 
CDC  6400,  or  about  0.1  seconds  on  a CDC  7600.  The  number  of  time  steps  required 
for  a given  computation  is  somewhat  less  than  those  required  for  the  nonlinear 
computations  at  small  values  of  K,  and  comparable  at  larger  values  of  K. 
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Conclusion 

An  accurate  and  efficient  procedure  for  computing  time-linearized,  small 
perturbation,  low  frequency  transonic  flows,  including  shock  wave  motions, 
has  been  developed.  Shock  motions  must  be  included  as  their  amplitude  is 
proportional  to  that  of  the  motion  divided  by  the  reduced  frequency.  Both 
indicial  and  harmonic  responses  for  various  modes  of  motion  may  be  computed 
in  seconds  on  a CDC  7600. 

f 
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FIGURE  8.  NORMALIZED  MAXIMUM  SHOCK  EXCURSION  AND 
CIRCULATION  AS  A FUNCTION  OF  INVERSE 
REDUCED  FREQUENCY  FOR  AN  NACA  64A006 
AIRFOIL  WITH  OSCILLATING  QUARTER-CHORD 
FLAP.  M_  - 0.375. 


